Software simulation tools allow physical systems to be studied without actually building the physical system. As an example, referring to FIG. 7, a two-dimensional physical system 204 includes a bead 200 moving on a vertically positioned wire 202 having a circular shape. System 204 can be described using equations:m{umlaut over (x)}=Fx,  (Equ. 1)mÿ=Fy−mg,  (Equ. 2)x2+y2=r2,  (Equ. 3)where x and y represent the coordinates of the bead, {umlaut over (x)} and ÿ represent the second derivative of x and y with respect to time t, respectively, m represents the mass of the bead, r represents the radius of the circle, and g represents the gravitational constant. Fx and Fy represent the forces exerted on bead 200 in the x and y directions, respectively. The motion of bead 200 on surface 202 can be determined by solving the equation:
                                                        [                                                                    m                                                        0                                                                                        0                                                        m                                                              ]                        ⁡                          [                                                                                          x                      ¨                                                                                                                                  y                      ¨                                                                                  ]                                =                                    [                                                                                          F                      x                                                                                                                                  F                      y                                                                                  ]                        +                                          2                ⁡                                  [                                                                                    x                                                                                                            y                                                                              ]                                            ⁢              λ                                      ,                            (                  Equ          .                                          ⁢          4                )            where λ is the Lagrangian multiplier. By differentiating Equ. 3 twice, we obtainx{umlaut over (x)}+yÿ+{dot over (x)}2+{dot over (y)}2=0.  (Equ. 5)Substituting Equs 1 and 2 into Equ. 5, we obtain
                              x          m                ⁢                  (                                    F              x                        +                          2              ⁢              x              ⁢                                                          ⁢              λ                                )                    +                        y          m                ⁢                  (                                    F              y                        -                          m              ⁢                                                          ⁢              g                        +                          2              ⁢              y              ⁢                                                          ⁢              λ                                )                    +                        x          .                2            +                        y          .                2              =    0    ,from which we obtain
                    λ        =                              -                                          m                ⁡                                  (                                                                                    x                        .                                            2                                        +                                                                  y                        .                                            2                                                        )                                                            2                ⁢                                  r                  2                                                              -                                    xF              x                                      2              ⁢                              r                2                                              -                                                    y                ⁡                                  (                                                            F                      y                                        -                                          m                      ⁢                                                                                          ⁢                      g                                                        )                                                            2                ⁢                                  r                  2                                                      .                                              (                  Equ          .                                          ⁢          6                )            Using Equs. 4 and 6, we can determine the positions of bead 200 over time using a discrete integration process (i.e., the position of the bead at time tn+1 is determined from its position at time tn by integrating differential terms).
Equs. 1-3 are said to be a set (or system) of differential algebraic equations (DAEs). In general, a system of DAEs can be written as:F(x′,x,y,t)=0,  (Equ. 7)G(x,y,t)=0,  (Equ. 8)where Equ. 7 represents one or more differential equations, Equ. 8 represents one or more algebraic constraints on the variables, and the Jacobian of F with respect to x′ (denoted by ∂F/∂x′=Fx′) is non-singular. By a change of variables, Equs. 7 and 8 can be rewritten into the form:F(t,u(t),u′(t))=0,  (Equ. 9)where the Jacobian of F with respect to u′ (Fu′) is singular.